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Hard ellipses: Equation of state, structure and self-diffusion 



<N 

o 

<N 
O 

Q 

00 
<N 



O 

s 

c 
o 
o 



Wen-Sheng Xu, Yan-Wei Li, Zhao-Yan Sun, and Li-Jia An 

State Key Laboratory of Polymer Physics and Chemistry, Changchun Institute of Applied Chemistry, 
Chinese Academy of Sciences, Changchun 130022, People's Republic of China 

(Dated: January 1, 2013) 

Despite their fundamental and practical interest, the physical properties of hard ellipses remain 
largely unknown. In this paper, we present an event-driven molecular dynamics study for hard 
ellipses and assess the effects of aspect ratio and area fraction on their physical properties. For 
state points in the plane of aspect ratio (1 < k < 9) and area fraction (0.01 < cj) < 0.8), we 
identify three different phases, including isotropic, plastic and nematic states. We analyze in detail 
the thermodynamic, structural and self- diffusive properties in the formed various phases of hard 
ellipses. The equation of state (EOS) is shown for a wide range of aspect ratios and is compared 
with the scaled particle theory (SPT) for the isotropic states. We find that SPT provides a good 
description of the EOS for the isotropic phase of hard ellipses. At large fixed 0, the reduced pressure 
p increases with k in both the isotropic and the plastic phases, and interestingly, its dependence 
on k is rather weak in the nematic phase. We rationalize the thermodynamics of hard ellipses 
in terms of particle motions. The static structures of hard ellipses are then investigated both 
posit ionally and orient at ionally in the different phases. The plastic crystal is shown to form for 
aspect ratios up to k = 1.4, while appearance of the stable nematic phase starts approximately at 
k = 3. We quantitatively determine the locations of the isotropic-plastic (I-P) transition and the 
isotropic-nematic (TN) transition by analyzing the bond-orientation correlations and the angular 
correlations, respectively. As expected, the I-P transition point is found to increase with k, while 
a larger k leads to a smaller area fraction where the I-N transition takes place. Moreover, our 
simulations strongly support that the two-dimensional nematic phase in hard ellipses has only quasi- 
long-range orientational order. The self- diffusion of hard ellipses is further explored and connections 
are revealed between the structure and the self- diffusion. We discuss the relevance of our results 
to the glass transition in hard ellipses. Finally, the results of the isodiffusivity lines are evaluated 
for hard ellipses and we discuss the effect of spatial dimension on the diffusive dynamics of hard 
ellipsoidal particles. 

PACS numbers: 61.20.Ja, 61.20.Lc, 61.30.Gd, 64.70.mf, 66.10.cg 



I. INTRODUCTION 



Hard-particle systems in general offer a great opportu- 
nity to understand the mechanism of a variety of physical 
phenomena, including crystal nucleation [HQ, glass tran- 
sition 0-0] and jamming (HQ. The spherical particles, 
such as hard spheres and hard disks, have been studied 
intensively during the past decades. In particular, for a 
system of hard spheres, an entropy-driven liquid-crystal 
transition occurs at sufficiently high density [loj and a 
hard-sphere glass can be also formed above the glass tran- 
sition point if crystallization is avoided [ll| . On the other 
hand, non-spherical particles play a fundamental role in 
the physics of molecular liquids and attract growing inter- 
est in recent years because of their usefulness in forming 
functional structures. From the theoretical point of view, 
it is interesting to include the rotational degrees of free- 
dom for the particles and investigate the physical behav- 
ior of a system composed of anisotropic particles. Such 
study uncovers new types of behavior and holds potential 
for providing new insights that cannot be obtained from 
studies of spherical objects [l2l-[l8j. From the practical 
point of view, the constituent particles often have non- 
spherical shapes in real materials such as laponite [l9| 
and the building blocks with anisotropic shape can pro- 
vide a powerful candidate for the assembly of particular 



targeted structures [20lJ2l|. Thus, exploration of systems 
composed of anisotropic particles is also valuable for ma- 
terials design. In this paper, we present an event-driven 
molecular dynamics (EDMD) study of hard ellipses. A 
system of hard ellipses can be regraded as one of the sim- 
plest models of two-dimensional (2D) anisotropic parti- 
cles and its phase behavior is very rich. Yet, the ther- 
modynamic, structural and dynamical properties of hard 
ellipses remain largely unknown. 

The phase transitions, in particular the isotropic- 
nematic (I-N) transition, of hard ellipses have been 
studied by computer simulations and various theoreti- 
cal approaches. In Vieillard-Baron's pioneering canon- 
ical Monte Carlo (MC) simulations [22|, three different 
phases, including isotropic liquid, nematic liquid crystal 
and solid, have been identified in a hard-ellipse system 
with an aspect ratio of k = a/b = 6, where a and b de- 
note semi-major and semi-minor axes. Vieillard-Baron 
observed that both I-N and nematic-solid (N-S) transi- 
tions are of first order, as evidenced by a small disconti- 
nuity shown in the equation of state (EOS, i.e., density- 
dependence of the reduced pressure p). However, com- 
puter simulations show that the I-N transition in a 2D 
fluid of hard needles [23[ and long hard rods Q , which 
in principle share the similar physics with elongated hard 
ellipses, is continuous rather than first-order. More- 
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over, the constant-pressure MC simulations of Cuesta 
and Frenkel [25| indicate that while there is no stable 
nematic phase for k = 2, the I-N transition is of first or- 
der for k = 4 and continuous via disclination unbinding 
for k = 6. Hence, it appears that the nature of the I-N 
transition depends on the aspect ratio and that the phase 
diagram possesses a tricritical point at which the transi- 
tion changes from first-order to continuous. Thus, on the 
simulation side, a definite conclusion has not be reached 
on the phase transitions of hard ellipses. This issue is also 
contentious on the theoretical side. Ward and Lado (26^ 
found based on the Percus-Yevick (PY) integral equation 
theory that the I-N transition does not take place in a 
fluid of hard ellipses although orientational ordering is al- 
lowed. While, the density-functional theory (DFT) [27| 
and the effective-liquid approach [28[ indicate that the I- 
N transition is continuous for hard ellipses with arbitrary 
aspect ratios. The other forms of DFT [29|, [30| and the 
scaled particle theory (SPT) (3l| predict the existence of 
the I-N transition in hard ellipses, but its location differs 
from one theory to another. Therefore, a clear under- 
standing of the phase transitions in hard ellipses is still 
lacking. Here, it should be also mentioned that the phase 
behavior of binary hard ellipses was recently investigated 
by a scaled particle DFT 132] and a weighted DFT [33|. 

Another important issue concerns the orientational 
order in the 2D nematic liquid crystals. Straley has 
shown [34[ that true long-range order (LRO) cannot exist 
in a 2D nematic phase if the particles interact via a sepa- 
rable potential, but the existence of true LRO cannot be 
excluded if the potential is not separable into positional 
and orientational parts. In fact, it has been demonstrated 
by computer simulations [23|, |24| that a general property 
of a 2D nematic phase is the lack of true LRO, although 
earlier simulations [35| also show that the system with 
the nonseparable potential exhibits true LRO. For this 
reason, the systems that lack true LRO are usually re- 
ferred to as having quasi-LRO. In the elastic continuum 
theory [36| , quasi-LRO is expected to appear in a 2D ne- 
matic phase if the free energy associated with collective 
fluctuations in the particle orientations can be expressed 
as 



(1) 



where 0(r) characterizes the orientation at position r 
with respect to a fixed axis and K is the 2D Frank's 
elastic constant. Based on Eq. (1), one can derive that 
the amplitude of the orientational fluctuations diverges 
with the system size in a logarithmic manner, 
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where < • • • > denotes the ensemble average, ks the 
Boltzmann's constant, T the absolute temperature and 
TV the particle number. As a consequence, both the 2D 
nematic order parameter q and the angular correlation 



function g^i decay algebraically, i.e., 

q =< cos(2#) >~ N~ kBT ^ K , 
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where I is a positive integer. Eqs. (3) and (4) imply that 
a 2D nematic phase with quasi-LRO can be characterized 
by a vanishing orientational order parameter in the ther- 
modynamic limit and a power-law decay of the angular 
correlation function. If the transition from a 2D nematic 
phase with quasi-LRO to an isotropic phase proceeds 
via the Kosterlitz-Thouless (KT) disclination unbinding 
mechanism [37j, it should occur at a universal value of 
the renormalized Frank's constant, 7rK c /8kBT = 1 (23| . 
As pointed out in Ref. [23] , a different mechanism is also 
possible, but there is no stable nematic phase at values 
of the Frank's constant below K c . Thus, one can quanti- 
tatively identify the location of the I-N transition if the 
nematic phase does form through the KT disclination un- 
binding mechanism. We note that the algebraic decay of 
the angular correlation function has been evidence d by 
simulations of 2D hard needles [23|, 2D hard rods [24[, 
hard ellipses [38], and recently by experiments of quasi- 
2D suspensions of hard ellipsoids [39| . 

In contrast to the efforts made on understanding the 
phase behavior and the structure of hard ellipses, little 
attention has been paid on their dynamics despite its fun- 
damental and technological importance. The dynamical 
properties of hard needles have been studied by molecular 
dynamics (MD) simulations about three decades ago (iof . 
The diffusive process of a single hard ellipsoid confined 
to two dimensions was measured just several years ago 
in experiments [4l|, |42|. In particular, the diffusion of 
quasi-2D suspensions of hard ellipsoids with k « 9 has 
been explored very recently by video-microscopy experi- 
ments [39| and the enhancement of the translational dif- 
fusion with respect to the rotational diffusion has been ra- 
tionalized in terms of the formation of unstable nematic- 
like regions with an average lifetime that exceeds the 
characteristic time of diffusion in a recent MC study [38| . 
In spite of the aforementioned progress, our knowledge of 
the dynamics of hard ellipses is far from complete. To our 
knowledge, there are no available MD results concerning 
the self-diffusion of hard ellipses so far. 

In this work, we analyze in detail the thermodynamic, 
structural and self- diffusive properties in the various 
phases of hard ellipses, with the aim of examining the 
effects of aspect ratio and area fraction and exploring 
connections between different properties. To this end, 
we perform EDMD simulations for hard ellipses over a 
wide range of aspect ratios and area fractions, covering 
state points including the isotropic, plastic and nematic 
phases. We present the EOS for a series of k values and 
compare the simulation results with a recent prediction 
of the SPT in the isotropic phase. We find that SPT 
describes well the isotropic branch of the EOS in hard 
ellipses. Dependence of k on the reduced pressure p is 
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explored in the various phases. At large fixed area frac- 
tions, p increases with k in both the isotropic and the 
plastic phases, and interestingly, its dependence on k is 
rather weak in the nematic phase. We rationalize the 
thermodynamics of hard ellipses in terms of particle mo- 
tions, which show different features in different phases. 
The static structures, both positionally and orientation- 
ally, are then investigated in the various phases. Our 
simulations show that the plastic crystal forms for as- 
pect ratios up to k = 1.4, while appearance of the stable 
nematic phase starts approximately at k = 3 for the stud- 
ied density range. The locations of the isotropic-plastic 
(I-P) transition and the I-N transition are determined 
by analyzing the bond-orientation correlation functions 
and the angular correlation functions, respectively. As 
expected, the I-P transition point is found to increase 
with /c, while a larger k leads to a smaller area fraction 
where the I-N transition takes place. In addition, our 
results strongly support that the 2D nematic phase in 
hard ellipses has only quasi-LRO. We further explore the 
self-diffusion of hard ellipses. We find that a phase tran- 
sition in the translational degrees of freedom is indeed 
reflected by an enhancement in the rotational diffusion 
and vice versa. Hence, a clear connection is revealed 
between the structure and the self-diffusion. We discuss 
the relevance of our results to the glass transition in hard 
ellipses, which is a subject of recent experimental inves- 
tigation [HI, EH- We finally provide the results of the 
isodiffusivity lines for hard ellipses and discuss the effect 
of spatial dimension on the diffusive dynamics of hard 
ellipsoidal particles by comparing results of hard ellipses 
and uniaxial hard ellipsoids [l3| . 

The paper is organized as follows. In Sec. II, we specify 
the details of the simulation and the model system used in 
this work. Thermodynamic, structural and self-diffusive 
properties of hard ellipses are presented and discussed 
in Sec. III. A summary of our findings is given in Sec. 
V. In the Appendix, we add the main ingredients for 
implementing an EDMD simulation for hard ellipses and 
a brief discussion on the effect of the system size. 



II. SIMULATION DETAILS 

We perform extensive EDMD simulations of hard el- 
lipses in a square box under periodic boundary condi- 
tions. The main ingredients for implementing an EDMD 
simulation for hard ellipses are described in the Appendix 
A. In a simulation, we fix the particle number N and 
the area of the simulation box V = L 2 with L the box 
dimension. In order to explore the properties in differ- 
ent phases, we investigate two sets of aspect ratios, i.e, 
k = 1 — 2 with an interval of 0.1 and k = 3 — 9 with 
an interval of 1, and a wide range of area fractions rang- 
ing from (j) = Nnab/V = 0.01 to 0.8. At sufficiently 
high densities, a nematic phase (i.e., particles have their 
centers of mass at random but exhibit some long-range 
orientational order) is expected to form for large k val- 



ues while a plastic phase (i.e., particles are oriented at 
random but their centers of mass have some ordered ar- 
rangement) can form for small k values. All the particles 
have the same mass m and the same moment of iner- 
tia /. Both m and / are set to be unity for convenience. 
The temperature T is irrelevant for athermal systems and 
remains constant due to the conservation of the total ki- 
netic energy. We set ksT to be unity in this work. Note 
that T contains both translational and rotational parts 
in the case of hard ellipses and that the separate tem- 
perature in each part fluctuates during the simulation. 
L ength, time and pressure are expressed in units of 26, 
Y // 4m6 2 /ksT and ksT/Ab 2 . In the following, we report 
results for N = 500. However, as can be seen from Eq. 
(3), it is also useful to study the effect of the system size 
especially in the nematic phase. We thus also perform 
simulations of hard ellipses with N = 200 and 800 for 
selected state points, and a brief discussion on the effect 
of the system size will be shown in the Appendix B. 

In order to generate the starting configuration of a 
hard-ellipse system with specific k and 0, we used the 
Lubachesky-Stillinger (LS) compression algorithm 0, 
|45[, in which the particles collide elastically and expand 
uniformly at a certain growth rate. LS algorithm was ini- 
tially designed for hard spherical particles and recently it 
has been generalized to the case of hard ellipsoidal parti- 
cles [46|, |4tJ. In the compression process, we started with 
a low-density hard-ellipse fluid at <p = 0.005 and used a 
quite small growth rate of 10 -4 so that the created initial 
configuration already remained close to be in an equilib- 
rium or stable state. These initial states were then used 
as inputs for EDMD simulations. At each aspect ratio, 
10 3 time units at <p < 0.6 and up to 10 4 time units at 
(j) > 0.6 were used to equilibrate the system, and then 
another 2 x 10 4 time units for collecting data. Four in- 
dependent runs were performed for each state point in 
order to obtain reliable results and improve the statis- 
tics. Thus, the sample- averaged results will be shown in 
the next unless otherwise stated. It is worth mentioning 
that the slowest run takes over four weeks on a 2.6 GHz 
CPU. 

We note that in the solid phase (i.e., particles show 
long-range order both positionally and orientationally ) , 
it will be useful to allow the shape of the simulation box 
to change during the simulation, because the box shape 
is necessarily compatible with the equilibrium shape of 
the unit cell in order to form a solid phase. This has been 
done in previous NPT MC simulations [25| in order to 
study the transition from an isotropic or a nematic phase 
to a solid phase in hard ellipses. The box shape is fixed to 
be square in our simulations since our aim here is not to 
determine the location where the solid phase forms, but 
to focus on the physical properties of the formed phases 
in EDMD simulations. Meanwhile, we do not investi- 
gate the melting process of a perfect hard-ellipse solid 
(which is a stretched triangular lattice). For these rea- 
sons, we didn't observe the solid phase for all the state 
points studied in this work, please see Subsection III B 



for further discussion. 



III. RESULTS AND DISCUSSION 

In this section, we present and discuss our results of 
the thermodynamic, structural and self-diffusive proper- 
ties in hard ellipses. First of all, the EOS is shown for a 
wide range of aspect ratios and the reduced pressure (i.e., 
compressibility factor) in the isotropic phase is compared 
with a recent result of the SPT. We also explore the de- 
pendence of the reduced pressure on the aspect ratio in 
the different phases and discuss the influence of the par- 
ticle motions on the thermodynamics. In the second sub- 
section, we analyze in detail the positional and orienta- 
tional order of hard ellipses and compare our results with 
previous simulations when possible. Finally, we investi- 
gate the self-diffusion of hard ellipses and demonstrate 
connections between the structure and the self-diffusion. 
The relevance of our results to the glass transition in hard 
ellipses is briefly discussed. We also show the results of 
the isodiffusivity lines in both the translational and the 
rotational degrees of freedom and discuss the effect of 
spatial dimension on the diffusive dynamics of hard el- 
lipsoidal particles. 



A. Equation of state 

Knowledge of an accurate EOS for a model system is of 
fundamental importance and particularly useful for pro- 
viding information on what type of an underlying phase 
transition is. We present the results of EOS for hard el- 
lipses at varying k in Fig. 1. The EOS for a system of 
hard disks is also included in Fig. 1(a). Here, the reduced 
pressure p = PV/NksT is computed by the momentum 
exchange from particle collisions in EDMD. We checked 
that our results for k = 2, 4 and 6 are consistent with the 
compression runs in previous MC simulations [22|, |25| . 

Let us first focus on hard ellipses with small k val- 
ues. In this situation, the orientational order is not ap- 
preciable in the system but the positional order can de- 
velop upon compression, leading to the formation of plas- 
tic crystals at sufficiently high densities. The appear- 
ance of the positional order can significantly affect the 
thermodynamics of the system and it is indeed reflected 
by a discontinuity in the EOS, as shown in Fig. 1(a). 
The plastic crystal was already observed for a system of 
hard ellipses with k = 1.01 about four decades ago by 
Vieillard-Baron [22|. Our results demonstrate that the 
plastic crystal can be formed in hard ellipses with aspect 
ratios up to k = 1.4. The observation of the disconti- 
nuity in the EOS is compatible with the scenario of a 
first-order transition, but it should be viewed with some 
caution since the nature of the formation of a 2D crys- 
talline phase is unclear even in the case of hard disks; 
e.g., the famous Kosterlitz, Thouless, Halperin, Nelson 
and Young (KTHNY) scenario [13, El, EH declares that 




FIG. 1: p as a function of (j> for (a) hard disks and hard 
ellipses with small k values and (b) hard ellipses with large 
k values. The results are shifted up by 5 from the preceding 
small k for clarity. The dotted lines are simulation results 
and the solid lines are the results of the SPT [Eq. (5)]. The 
green dashed line in (a) is the SPT prediction for the equation 
of state of hard disks in the isotropic state, i.e., Eq. (10) 
in Ref. [53|]. The circles indicate the area fractions where 
the system starts to form stable plastic crystals in (a) and 
stable nematic crystals in (b), respectively (see the structural 
analysis in the next subsection). The semi- log plot in the 
inset of Fig. 1(b) highlights the change of dependence of p on 
(f) with increasing k. 



the melting of hard disks proceeds via two continuous 
phase transitions with an intermediate phase called hex- 
atic phase, while the recent large-scale computer sim- 
ulations [50] reveal that the liquid-hexatic transition is 
of first order. We will further discuss this issue in the 
next subsection by analyzing the structure of the sys- 
tem. Moreover, we find in Fig. 1(a) that the disconti- 
nuity shown in the EOS becomes less apparent as k gets 
larger and nearly undetectable within the studied den- 
sity range for k > 1.5. This can be easily understood 
since the particle's anisotropy tends to destroy the posi- 
tional order of the system, and thus will be in favor of the 
isotropic phase. For hard ellipses with sufficiently large 
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FIG. 2: p as a function of 7 at fixed densities. The solid lines 
are the results of Eq. (5). 



elongations, the orient at ional order is expected to appear 
and the positional order will be highly suppressed, result- 
ing in the formation of 2D nematic liquid crystals in hard 
ellipses. As already found in other 2D models [23|, l24j , 
the thermodynamic properties of the system are indeed 
not very sensitive to the occurrence of the orientational 
order. This is also confirmed by our results of the EOS 
for hard ellipses with k > 3, as shown in Fig. 1(b). For 
these k values, the system does have a nematic phase 
above some area fraction, but the corresponding change 
in the EOS is rather mild. Thus, we cannot unambigu- 
ously determine the nature of the I-N transition based 
on solely the thermodynamics. However, by a close in- 
spection of the results in the semi-log plot [see the inset 
of Fig. 1(b)], some difference in the EOS can be seen 
with varying k. It seems that some inflection points ex- 
hibit in the EOS for large aspect ratios, while the slope 
d\ogp/d(j) monotonically increases with <\> for small as- 
pect ratios, i.e., the ^-evolution of p differs qualitatively 
from relatively small to very large aspect ratio. This sug- 
gests that the nature of the I-N transition changes with 
aspect ratio. In fact, previous MC simulations imply the 
existence a tricritical point between k = 4 and 6 in the 
phase diagram of hard ellipses [25| . The I-N transition is 
first-order below the tricritical point and becomes contin- 
uous above it. Interestingly, we find that the appearance 
of the inflection points in the EOS starts within the same 
k range. Thus, our results appear to support the exis- 
tence of a tricritical point in the phase diagram of hard 
ellipses. However, since the location of such a point is 
very hard to determine even if it exists, we will focus on 
characterizing the properties of different phases in hard 
ellipses. Next, we compare the reduced pressure in the 
isotropic phase of hard ellipses with a recent result of the 
SPT. 

As mentioned in Section I, the EOS of hard ellipses 
can be also obtained by various theories. Hence, it is 
meaningful to test these theoretical predictions as our 



results cover a wide window of k and <p and contain data 
in the various phases. SPT in general provides a simple 
form for the EOS of a model system and its usefulness has 
been confirmed by the study of hard disks [5l|, [52| . Other 
existing theories also provide relevant results, but the 
solution of the EOS needs to be solved numerically (27L 
I29I43H . Thus, we concentrate here on a very recent result 
of the SPT [HI. On the basis of the SPT, Boulik derived 
an EOS for the isotropic fluid of hard ellipses [53|, which 
has the following form 
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where the non-circularity parameter 7 has been intro- 
duced. 7 is calculated according to the formula 7 = 
C 2 1 (4tt 2 ab) with C the perimeter of the ellipse. The 
value of C can be expressed accurately by the complete 
elliptic integral of the second kind, and a rather good ap- 
proximation can be obtained from the following formula: 
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(6) 



where v = (a — b)/(a+b). We use Eq. (6) to calculate 7 in 
this work and the results of the Eq. (5) are shown as solid 
lines in Fig. 1. As a reference, the area fractions where 
the system starts to form stable plastic crystals and sta- 
ble nematic crystals are also shown in Fig. 1. As can 
be seen in Fig. 1(a), the simulation results and the SPT 
prediction are almost indistinguishable in the isotropic 
phase for k < 2, indicating a fairly good description of 
Eq. (5) for the EOS of hard ellipses with small k values. 
Meanwhile, the predicted EOS deviates from the simula- 
tion results above some area fraction for k > 3, as shown 
in Fig. 1(b). As explained in Ref. [53| . the system is not 
isotropic but exhibits orientational order at high densi- 
ties for large k values. In Fig. 1(b), one can find that 
the area fraction where Eq. (5) starts to deviate from the 
simulations is much smaller than the I-N transition point 
[see circles in the main plot of Fig. 1(b)]. However, it 
should be emphasized that the short-range orientational 
order already appears in the system at smaller area frac- 
tions (see next subsection). Taking this into account, we 
find that Eq. (5) describes well the isotropic branch of 
the EOS for hard ellipses even with large k values. We 
can also use 7 as a variable to test Eq. (5). The results 
for representative area fractions are shown in Fig. 2. We 
find that the theoretical and simulation results agree well 
for low densities. At high densities, deviations appear for 
large 7 values because of the occurrence of short-range or 
some long-range orientational order, as evidenced by the 
non-monotonic dependence of p on 7. Therefore, we con- 
clude that SPT provides a good description of the EOS 
for the isotropic phase of hard ellipses. 

It is interesting to further explore dependence of the 
reduced pressure on the aspect ratio in the various phases 
of hard ellipses. To this end, we plot p as a function of 
k for (j> = 0.78, 0.79 and 0.8 in Fig. 3. For these den- 
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FIG. 3: /c- dependence of p at the highest densities studied. 
The dotted lines are guide to the eye. For these densities, the 
system has a plastic phase for k < 1.4, a nematic phase for 
k > 3 and an isotropic phase in between. 



sities, the system has a plastic phase for k < 1.4, a ne- 
matic phase for k > 3 and an isotropic phase in between 
(see next subsection). We find that p monotonically in- 
creases with k at fixed densities in both the isotropic 
and the plastic phases (i.e., k < 2 at <p > 0.78), sug- 
gesting that the thermodynamics of the system can be 
significantly influenced by the particle shape when the 
ellipses are oriented at random. Note that this conclu- 
sion is still valid for even larger k values (see the results 
for cj) < 0.2 in Fig. 2). Interestingly enough, in the ne- 
matic phase (i.e., k > 3 at <p > 0.78), when the particles 
show some long-range orientational order but have their 
centers of mass at random, p depends on k very weakly. 
Although not pointed out before, the same conclusion 
can be drawn from previous NPT MC simulations [25^ . 
which indicate that the area fractions are almost iden- 
tical at the same pressure in the nematic phase of hard 
ellipse at least for k — 4 and 6. To understand the above 
findings, we recall that the reduced pressure is calculated 
from the momentum exchange during collisions of two el- 
lipses, which indeed contains translational and rotational 
parts (although they are coupled with each other). It is 
expected that the contribution from the rotational mo- 
tions becomes more apparent as the ellipses get more 
elongated. This is the reason why p increases with k in 
the plastic and the isotropic phases, where the particles 
can freely move both translationally and rotationally (see 
Fig. 14). In the nematic phase, however, particles on av- 
erage line up preferentially along a common direction. 
In this case, the rotational dynamics becomes extremely 
slow, but the motions of the particles have little hin- 
drance in the translational degrees of freedom (see Fig. 
14) and thus the translational part becomes dominant 
in determining the reduced pressure of the system. We 
then speculate that the translational part must be very 
similar for all aspect ratios at the same density in the 




FIG. 4: g(r) of hard ellipses at <j> = 0.8 for k > 1.5. The 
results are shifted up by 0.5 from the preceding small k for 
clarity, As there are two basic length scales in a system of 
hard ellipses, g(r) also peaks at r = 2a, as indicated by the 
arrows for k — 9. 



nematic phase so that p has only weak dependence on k. 
Therefore, the thermodynamics of hard ellipses could be 
explained in terms of the particle motions. 



B. Static structure 

We now turn to the structural properties of hard el- 
lipses. As a system of hard ellipses with sufficiently large 
aspect ratios has a transition from the isotropic liquid to 
the 2D liquid crystal, much attention has been paid on 
the structure of the nematic phase in hard ellipses, i.e., on 
characterizing the orientational order of elongated hard 
ellipses. Instead of focusing solely on the orientational 
correlations, we explore both positional and orientational 
order in the various phases of hard ellipses in this sub- 
section. 

We first consider the pair correlation function of hard 
ellipses, which is defined as 



](r) 



27rrArN(N - 1) 



<^(r-M)>, (7) 



where r is the distance between the centers of mass of 
the particles. In general, characteristic peaks will appear 
in g(r) at large distances if the system possesses some 
long-range positional order. This is indeed observed in 
systems of hard ellipses with k < 1.4 at high densities 
(data not shown), where the centers of mass of the par- 
ticles form a triangular (or hexagonal) lattice. In hard 
ellipses with larger aspect ratios, an ordered structure in 
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FIG. 5: gQ(r)/g(r) of hard ellipses with k = 1.2 in the vicin- 
ity of the isotropic-plastic transition. The solid lines are the 
results of the OZ fittings and the dash dotted lines are the 
results of the power-law fittings. The green dashed line in- 
dicates gQ(r)/g(r) ~ r _1//4 . The results are similar for other 
aspect ratios with k < 1.4. 



principle can also form translationally although it needs 
not to be a perfect triangular one. We present the results 
of g(r) at <p = 0.8 for k > 1.5 in Fig. 4. Note that there 
are two basic length scales (i.e., the semi- major axis a 
and the semi-minor axis b) in a system of hard ellipses 
so that g(r) also peaks at the distance of 2a. Beyond 
the distance of 2a, we find from Fig. 4 that g(r) rapidly 
decays to the value of 1, indicating the absence of any 
long-range translational order even at the largest density 
studied. The reason, why we didn't observe the forma- 
tion of a translationally ordered structure for large aspect 
ratios, has been explained in Section II. Moreover, it is 
clear that several peaks are also shown in g(r) at between 
2b and 2a, suggesting the existence of some short-range 
translational order. As we shall see later, this is indeed 
due to the appearance of the orient ational order. 

The pair correlation function only provides a qualita- 
tive way to see whether a system possesses long-range 
positional order or not. Meanwhile, we can study the 
bond-orientation correlation functions [54|, in order to 
quantitatively identify the location of the transition from 
the isotropic liquid to the hexagonal crystal. We first de- 
fine the sixfold bond-orientation order parameter 

rij 

4 = ~J2 eMiWl), (8) 

Tij 
J m—1 

where i = \/— 1, rij is the number of the neighbors of 
the jth particle and 6 J m is the angle between the vector 
(r m — Yj) and the x axis. Here, we define two ellipses to 
be neighbors if they overlap when uniformly expanding 
their sizes by 1.4. Then, the spatial correlation of tp J 6 can 



be calculated by 

*M = 2 ,rArN(N-l) < 5> " > ' 

(9) 

The quasi-long-range bond-orientation order is evidenced 
by an algebraic decay of gsir) / 'g(r), i.e., 

9%{r)/g{r)^r-^. (10) 

In particular, according to the KTHNY theory [48|, t]q 
has a value of 1/4 at the boundary between the liquid 
phase and the hexatic phase. We thus use the criterion 
gs(r)/g(r) ~ r -1 / 4 to separate the plastic states from 
the isotropic states for hard ellipses with k < 1.4. In 
the liquid phase, the system only exhibits short-range 
bond-orientation order and g§(r)/g(r) should decay ex- 
ponentially. However, for high densities and especially 
in vicinity of the formation of the hexatic phase, we find 
that the bond-orientation correlation function can indeed 
be better described by the Ornstein-Zernike (OZ) func- 
tion 

^6(r)/^(r)^r- 1 / 2 exp(-r/e 6 ), (11) 

with ^6 the static correlation length. In fact, it is widely 
observed that the OZ function can characterize well the 
behavior of g§{r)/g{r) in 2D supercooled liquids p55l-[57|. 
We confirm the above claims by presenting the results of 
96(r)/g(r) for a system of hard ellipses with k = 1.2 
in Fig. 5. As can be seen, g§{r)/g{r) decays slowly 
with increasing (j). Below (j) ~ 0.72, the envelopes of 
gQ{r)/g(r) can be well fitted by the OZ function. With 
further increasing (/>, a power-law decay with an exponent 
of t)q < 1/4 holds for the behavior of gsir) / g{r), indicat- 
ing the emergence of quasi-long-range bond-orientation 
order. We thus estimate the location of the isotropic- 
plastic transition to be <p p = 0.72 for k = 1.2. Note that 
(j) p obtained in this way is actually an upper bound to 
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FIG. 6: cj) p as a function of k. The solid line is a guide to the 
eye. 
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FIG. 7: P2 as a function of Odir for hard ellipses with k — 6 at 
varying 0. Note that the result is shown for a single sample 
because P2(0dir) cannot be averaged among different samples 
since the nematic direction differs from one sample to another. 

the exact value. The values of (j) p for other aspect ratios 
are shown in Fig. 6. and we confirmed that the obtained 
value for hard disks (i.e., k = 1) is agreement with other 
studies within the simulation accuracy [Ho|. The above 
analysis indicates that hard ellipses for k < 1.4 share 
the similar behavior with hard disks on the structure as 
well as the thermodynamics (see discussion in the preced- 
ing subsection). Thus, the nature of the isotropic-plastic 
transition of hard ellipses should be same as that of hard 
disks. 

Turing to the orient ational order of hard ellipses, we 
first focus on the nematic order parameter as a function 
of cj) and k. The 2D nematic order parameter is defined 
as 

1 N 

P2 = jj < ^2cos 2 (^ - dir ) > -1, (12) 

j=l 

where 0j is the angle characterizing the orientation of the 
jth ellipse with respect to the x axis and Odir the orien- 
tation of the nematic director. Note that this definition 
is equivalent to Eq. (3). As the nematic director is not 
known a priori, P2 and Odir are usually determined by 
finding the eigenvalues and the eigenvectors of a tensor 
order parameter I23l 12511. As pointed out and observed 
in previous work [23|, [25j , the nematic order parameter 
obtained from the above method has a clearly nonzero 
value even in the isotropic phase due to the finite system 
size used in simulations. Instead of using the tensor or- 
der parameter, we set Odir as a variable and calculate P2 
as a function of Odir in this work. We show the results of 
P2(6dir) at varying <p for k = 6 in Fig. 7 as an illustra- 
tion. As can be seen, P2(0dir) shows extrema, and the 
absolute values of these extrema are the same (we de- 
note it by p™ ax ) because of the periodicity of the cosine 
function in the definition of P2. We then use P™ ax to 
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FIG. 8: Upper: p™ ax as a function of (j) at varying k. Lower: 
contour plot of p^ ax in the plane of (f> and k. 



quantify the nematic order of hard ellipses. In Fig. 8, 
pmax p reseri ted as a function of (j) at various aspect 
ratios and also shown as a contour plot in the plane of 
(j) and k. We find from Fig. 8 that P™ ax is nearly zero 
at low densities and has clearly nonzero values at high 
densities for k > 3. Thus, it appears that P™ ax is more 
sensitive to the onset of the orientational order than q. 
For k < 2, a clearly nonzero but still very small value 
of P™ ax sets in above <fi w 0.7, implying the absence of 
any long-range orientational order. The quasi- long-range 
orientational order is expected to appear for systems of 
hard ellipses with larger aspect ratios, as confirmed by 
the results of P™ ax for k > 3. Moreover, the area frac- 
tion where the system starts to orientate becomes small 
as k increases, revealing the fact that a system of more 
elongated hard ellipses has a stronger tendency to form 
nematic liquid crystals. 

It is interesting to ask whether P™ ax is dependent on 
the system size in the nematic phase. Unlike the order 
parameter q computed from the tensor order parameter 
(dependence of the system size on q has been confirmed 
by MC studies [23-25]), we find that the system size has 
little influence on P™ ax in both the isotropic and the ne- 
matic phases (see Appendix B). Hence, it is not proper 
to estimate the point of the I-N transition from the de- 



9 




FIG. 9: #2(7*) for hard ellipses with k < 2 at 
green dashed line indicates #2(7) ~ r -1 ^ 4 . 



0.8. The 



pendence of the system size on P™ ax . We can quanti- 
tatively determine the location of the I-N transition <j) n 
from the angular correlation function defined by Eq. (4). 
We mainly concentrate on the case of I = 1, i.e., 



g 2 (r) =< cos(2[6>(r) -0(0)]) > 



no r 



(13) 



where 7/2 = IksT ji\K and the average is performed over 
all pairs with the distance r. As introduced in Section I, 
the stable nematic phase is then identified by a power- 
law decay with 772 < 1/4 for #2(7*)- The results of g^ir) 
are shown for k < 2 at = 0.8 in Fig. 9. Consistent 
with the results of P™ ax ? we find that gii?) decays much 
faster than giif) ~ r -1 / 4 even at the largest studied den- 
sity for these aspect ratios, indicating again there is no 
quasi-long-range orientational order in the system. The 
algebraic decay appears at high densities for k > 3, as 
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FIG. 10: g2(r) for a system of hard ellipses with k = 6 in the 
vicinity of the isotropic-nematic transition. The solid lines 
are the results of the power-law fittings. The green dashed 
line indicates g2(r) ~ r _1//4 . The results are similar for other 
aspect ratios with k > 3. 
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FIG. 11: cj)n as a function of k. The triangle and the squares 
are the MC results taken from Ref . [22|] and Ref . [25|] , respec- 
tively. The solid line is a guide to the eye. 



evidenced in Fig. 10 for a system of hard ellipses with 
k = 6. We then use the power law #2(7") ~ r ~ m to fit 
g2(f) at r > 2a and estimate the location of the I-N tran- 
sition (j) n to be the area fraction when the value of 772 is 
smaller than 1/4. Again, </> n obtained in this way is actu- 
ally an upper bound to the exact value. It should be also 
stressed that the I-N transition in a finite system tends 
to occur at a lower density than in an infinite system be- 
cause the Frank's constant obtained from simulations on 
small systems is larger than that in the infinite system 
size [23]. We show the results of (j) n as a function of k in 
Fig. 11. As can be expected, <j) n decreases with increas- 
ing k. Comparing with previous MC results, we find that 
(j) n for k = 6 obtained in this work is larger than that in 
Ref. [22] but close to the result in Ref. (25j |. The reason, 
why the I-N transition was observed at a smaller area 
fraction by Vieillard-Baron, has been hinted by Frenkel 
and Eppenga [23[. Similar to the case of k = 6, the value 




FIG. 12: g 2 (r) for hard ellipses with k > 3 at = 0.1 
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of <j) n for k = 4 in this work is also a little smaller than 
that of Ref. 12511 . The slight difference between the re- 
sults of Ref. [25j and ours is possibly due to the different 
particle number and different method used (N « 200 and 
NPT MC in Ref. [H versus N = 500 and NVT MD in 
this work). We also note a recent MC study [38|, where 
(j) n for k = 9 is roughly estimated to be 0.58, which is 
much larger than ours. On the other hand, as mentioned 
in Section I, the nature of the orientational order in the 
2D nematic phase is also a central theme in the study of 
liquid crystals. The algebraic decay of #2(7) already in- 
dicates the lack of the true long-range orientational order 
in the nematic phase of hard ellipses. We further confirm 
that the algebraic orientational order holds even for the 
largest studied density, as shown in Fig. 12, although 
#2(0 ma y decay very slowly in this case. Thus, our sim- 
ulations strongly support that the nematic phase in hard 
ellipses has quasi-LRO. 

Based on the above structural analysis, we present a 
phase diagram of hard ellipses and several representative 
snapshots in Fig. 13. Our simulations show three dif- 
ferent phases within the investigated state points. For 
k < 1.4, the system has a transition from a low-density 
isotropic liquid to a plastic crystal at high densities [Fig. 
13(d)]. With increasing the aspect ratio up to k = 2, 
we didn't observe any phase transition within the whole 
density range studied in this work [Fig. 13(b)]. For hard 
ellipses with k > 3, an isotropic phase [Fig. 13(c)] will 
transform into a nematic phase [Fig. 13(e)] at sufficiently 
high densities. Again, we confirm from the snapshot in 
Fig. 13(c) that some short-range orientational order ex- 
ists even in the isotropic phase of elongated hard ellipses. 
On the other hand, although the simulated phase dia- 
gram does not include any solid phase, we should stress 
again that hard ellipses do have the transition from an 
isotropic or a nematic phase to a solid phase, as revealed 
in previous MC simulations [22|, [25j. We even expect 
that a transition from a plastic crystal to a solid phase 
will occur for small k values at even larger densities than 
(j) = 0.8. The reason why such a transition was not ob- 
served in our work, as explained in Section II, is due to 
the fact that a solid is hard to spontaneously form from a 
liquid in NVT MD simulations. In fact, our knowledge is 
very limited so far about the equilibrium structure of the 
solid in hard ellipsoidal particles, although it can be read- 
ily shown that the densest packing of hard ellipses has the 
same density as the densest packing of hard disks [l5| . 
More elegant methods will be needed to clarify this issue. 
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FIG. 13: (a) Phase diagram of hard ellipses in the plane of k 
and <j). Snapshots of (b) a high-density isotropic phase with 
k = 1.5 and 4> = 0.8, (c) a low-density isotropic phase with 
k = 6 and 4> = 0.5, (d) a plastic phase with k = 1.2 and 
cj) — 0.8, (e) a nematic phase with k — 9 and <j) — 0.8. Note 
that there is visually no distinction for the orientation of an 
ellipse with and + tt so that particles with is shown as 
the same color as those with + tt in the snapshots. 



C. Self-diffusion 

Self-diffusion is of great importance in many real pro- 
cesses, but it is still poorly understood in hard ellipses. 
The diffusive property of hard ellipses with k = 9 has 
been studied recently by experiments [39| and MC simu- 
lations (38[. Yet, the effects of aspect ratio and area frac- 
tion on the self-diffusion of hard ellipses are unclear. In 



this subsection, we investigate in detail the self-diffusion 
of hard ellipses over a wide range of aspect ratios and 
area fractions, including the isotropic, plastic and ne- 
matic phases, and demonstrate the connections between 
the structure and the self-diffusion. We further discuss 
the relevance of our results to the glass transition in hard 
ellipses and provide the results of the isodiffusivity lines. 

The self-diffusion is measured here by the translational 
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FIG. 14: Time evolution of (a) translational and (b) rota- 
tional MSDs at cf) = 0.8 for three aspect ratios. The system 
has a plastic phase for k = 1.1, an isotropic phase for k = 2 
and a nematic phase for k = 6. 



FIG. 15: (a) Translational diffusion constant Dt and (b) ro- 
tational diffusion constant Dq as a function of <j> at varying 
k. 



and rotational mean squared displacements (MSD) 
1 N 

<r 2 (t)>=-<^|r J (t)-r J (0)| 2 >, (14) 



N 



< e\t)>=^<Y j \e J (t)-0M\ 2 >, (is) 



and the corresponding diffusion constants obtained from 
the Einstein relation 



D r = lim 



< r 2 (t) > 



t^oo At 



Do 



lim 

t— >-oo 



< e 2 (t) > 

2t 



(16) 



(17) 



We first focus on the behavior of the MSD in the differ- 
ent phases of hard ellipses. The results of < r 2 (t) > and 



< 9 2 {t) > are shown at = 0.8 for the plastic (k = 1.1), 
isotropic (k = 2) and nematic (k = 6) phases of hard 
ellipses in Fig. 14. As can be seen, both translational 
and rotational motions of particles are ballistic at short 
times, and crosses over into diffusive behavior at suffi- 
ciently long times in both the plastic and the isotropic 
phases. Moreover, an apparent subdiffusive behavior is 
seen in the intermediate time regime of the translational 
MSD due to the cage effect in the plastic phase [see the 
solid line in Fig. 14(a)]. By contrast, the rotational dif- 
fusion of particles becomes extremely slow in the nematic 
phase, while the translational MSD can still retain a dif- 
fusive regime at long times within the simulation time 
window (see the dotted lines in Fig. 14). We have seen 
in Subsection A that the above features shown in the par- 
ticle motions may be responsible for the thermodynamics 
in the different phases of hard ellipses. 

To better illustrate the effects of aspect ratio and area 
fraction on the self-diffusion of hard ellipses, we extract 
the translational and rotational diffusion constants from 
those of the corresponding MSDs retaining a diffusive 
regime. The results of Dt and Dq are shown in Fig. 
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FIG. 16: Dq as a function of Dt at varying k. The green 
dashed line indicates Dq — Dr- 



ib. At fixed area fractions, we observe that increasing k 
leads to a slowing down of the rotational diffusion, while 
Dt increases with aspect ratio. Moreover, the occur- 
rence of the phase transition seems not to affect these 
trends at least within the studied density range. Then, 
a natural question to ask is how self-diffusion responds 
to the phase transition. In a recent MD study j38|, it 
was found for a system of hard ellipses with k = 9 that 
if the diffusion constants are normalized by the values in 
the infinite dilution limit, then both diffusion constants 
still monotonically decrease with <\> but the normalized 
Dt will exceed the normalized Dq at the area fraction 
where the I-N transition takes place, i.e., the transla- 
tional diffusion is enhanced with respect to the rotational 
diffusion as a result of the appearance of the quasi-long- 
range orientational order. We didn't observe crossing of 
the two normalized curves for the same aspect ratio by 
performing the same analysis (data not shown). Instead, 
our simulations reveal that the self-diffusion responds to 
the phase transition in a rather straight manner, which 
can be detected already in each curve; i.e., a sudden in- 
crease of Dq is seen when a plastic crystal forms, while 
the emergence of a nematic phase leads to an enhance- 
ment of Dt although this observation seems to become 
less evident as k gets smaller. Hence, a phase transition 
in the translational degrees of freedom is reflected by a 
corresponding change in the rotational diffusion and vice 
versa. We thus uncover a clear connection between the 
structure and the self-diffusion in hard ellipses. It is in- 
teresting to further explore the effect of the particle's 
anisotropy on the self-diffusion in each degrees of free- 
dom by plotting Dq as a function of Dt at varying k. As 
shown in Fig. 16, Dq is larger than Dt within the full 
density range for k < 2 and the particles diffuse slower 
in the rotational degrees of freedom than in the trans- 
lational degrees of freedom for k > 3. In addition, the 
curves show an upturn at high densities for k < 2, while 
they turn down for k > 3. This means that the trans- 
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FIG. 17: Isodiffusivity lines in the plane of area fraction and 
aspect ratio. The solid lines are isodiffusivity lines from trans- 
lational diffusion coefficients Dt and the dashed lines are isod- 
iffusivity lines from rotational diffusion coefficients Dq. The 
green squares and the red circles indicate the locations of the 
isotropic-plastic transition and the isotropic-nematic transi- 
tion, respectively. 

lational mobility of hard ellipses will be similar with the 
rotational mobility at some aspect ratio between 2 and 
3. If a glass is allowed to form in hard ellipses, we then 
expect that the rotational glass transition sets in at a 
lower density than the translational glass transition for 
k > 3, which was indeed observed recently in experiments 
of quasi-2D hard ellipsoids with k ~ 6 [18]. Instead, the 
translational glass transition density will be smaller than 
the rotational one for k < 2. This immediately suggests 
that there is a point in the phase diagram of glassy hard 
ellipses, where the translational and the rotational glass 
transition lines will intersect. Therefore, glasses in hard 
ellipses can be classified into three categories: a plastic 
glass in which dynamic arrest occurs in the translational 
degrees of freedom but the system remains a liquid in the 
rotational degrees of freedom, a liquid glass in which the 
system becomes glass in the rotational degrees of free- 
dom but the translational motions of particles remains 
ergodic and a glass in which motions of the particles are 
dynamically arrested in both the translational and the 
rotational degrees of freedom. 

Another relevant property is the isodiffusivity of hard 
ellipses. As there are related results for hard ellipsoids of 
revolution [13|, we can evaluate such isodiffusivity lines 
for hard ellipses by proper interpolation and thus explore 
the qualitative effect of spatial dimension. We present 
the results of the isodiffusivity lines in Fig. 17. We 
find that an almost perpendicular crossing of the trans- 
lational and the rotational isodiffusivity lines is shown 
in the plane of (j) and k and that the rotational isodif- 
fusivity lines of hard ellipses reproduce qualitatively the 
shape of the I-N transition line at least for sufficiently 
small Dq values, findings that are similar with the results 
in hard ellipsoids of revolution [l3|. The latter observa- 
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tion indicates that the rotational motions of the ellipses 
is mostly controlled by the particle's anisotropy. On the 
other hand, the translational isodiffusivity lines in uni- 
axial hard ellipsoids mimic the swallowlike shape of the 
coexistence line between the isotropic liquid and the crys- 
talline phase [13| , suggesting a non-monotonic change of 
Dt with k at fixed density (see Fig. 2 in Ref. [l3j). 
Note that the swallowlike shape appears because there is 
a distinction between prolate (rodlike) and oblate (disk- 
like) particles for hard ellipsoids, while such a distinction 
does not exist in hard ellipses since there are only two 
symmetry axes. Our results, already apparent from Fig. 
15, show that Dt monotonically increases with k at fixed 
(/>, which remains valid even when a phase transition in- 
tervenes. Thus, it appears that the spatial dimension 
does influence the translational dynamics. However, it 
should be noted that our result does not rule out a non- 
monotonic dependence of Dt on k at even larger fixed 
densities than that presented in Fig. 17. In fact, our 
preliminary study for polydisperse hard ellipses (which 
prevent the formation of crystal phases) does show that 
the translational glassy dynamics depends on aspect ra- 
tio non- monotonically at large fixed densities. Therefore, 
a common mechanism will be responsible for the trans- 
lational glassy dynamics of hard ellipsoidal particles al- 
though different origins may exist at low densities. 



IV. CONCLUSIONS 

In summary, we have presented a numerical study of 
hard ellipses over a wide range of aspect ratios and area 
fractions, covering state points including the isotropic, 
plastic and nematic phases. We have provided the re- 
sults of the EOS for a series of k values and compared 
the simulations with a recent prediction of the SPT in 
the isotropic phase. We find that SPT describes well 
the isotropic branch of the EOS in hard ellipses. Depen- 
dence of k on the reduced pressure p has been explored 
in the various phases. At large fixed area fractions, p in- 
creases with k in both the isotropic and the plastic phase, 
and interestingly, its dependence on k is rather weak in 
the nematic phase. We rationalize the thermodynam- 
ics of hard ellipses in terms of particle motions, which 
show different features in different phases. For the static 
structures, our simulations show that the plastic crys- 
tal forms for aspect ratios up to k = 1.4, while appear- 
ance of the stable nematic phase starts approximately 
at k = 3. The locations of the I-P transition and the 
I-N transition have been determined respectively by ana- 
lyzing the bond-orientation correlation functions and the 
angular correlation functions and compared with previ- 
ous simulations. We demonstrate that the I-P transition 
point increases with k as the particle's anisotropy tends 
to destroy the positional order, while a larger k leads to a 
smaller area fraction where the I-N transition takes place, 
revealing the fact that a system of more elongated hard 
ellipses has a stronger tendency to form nematic liquid 



crystals. In addition, our results strongly support that 
the 2D nematic phase in hard ellipses has only quasi- 
LRO, as evidenced by the algebraic decay of the angular 
correlation function in the nematic phase. We have inves- 
tigated the self-diffusion of hard ellipses and found that 
a phase transition in the translational degrees of freedom 
is indeed reflected by a corresponding change in the rota- 
tional diffusion and vice versa. Hence, we reveal a clear 
connection between the structure and the self-diffusion. 
We discuss the relevance of our results to the glass transi- 
tion in hard ellipses and infer that there is an intersection 
point between the translational and the rotational glass 
transition lines in the phase diagram of glassy hard el- 
lipses. Finally, we have evaluated the isodiffusivity lines 
for hard ellipses and discussed the effect of spatial di- 
mension on the translational dynamics by comparing the 
diffusive dynamics of hard ellipses with that of hard el- 
lipsoids. Our results are also valuable for understanding 
the structure and the diffusion of anisotropic molecules 
at membranes and interfaces. 
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Appendix A: Event-driven molecular dynamics for 
hard ellipses 

We describe the main ingredients for implementing an 
EDMD simulation of hard ellipses in this section. In 
order to perform a MD or a MC simulation for a hard- 
particle system, a key task is to detect overlap between 
two particles. While this is a trivial thing for spherical 
particles, it is highly non-trivial for non-spherical parti- 
cles. In the case of hard ellipses, Vieillard-Baron [22j first 
introduced an overlap criterion about four decades ago. 
Later on, Perram and Wertheim [58| derived a contact 
function. Both criteria can be easily applied in a numer- 
ical simulation of hard ellipses. However, it appears that 
the Perram- Wertheim (PW) approach is more convenient 
both computationally and theoretically. Moreover, in ad- 
dition to an overlap criterion, we will need the first-order 
time derivative of the distance of closest approach be- 
tween two moving ellipses in order to determine the con- 
tact time, which is a central quantity in an EDMD sim- 
ulation. Indeed, Donev et al [47} have derived the time 
derivatives of the PW contact function and developed an 
EDMD algorithm for hard ellipsoidal particles. We thus 
follow their results and adopt the PW approach in our 
study. However, we note that other approaches also exist 
for detecting overlap between two ellipses (591461) . 

Let us first describe the orientation of a rigid body. In 
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two dimensions, it needs only one angle to represent the 
orientation of the rigid body as there is only one rational 
degree of freedom. Let 6 denote the angle between the 
semi-major axis of the ellipse and the x axis, then the 
rotation matrix reads 



Q = 



cos sin 
— sin 9 cos 6 



(Al) 



Taking its shape into account, an ellipse can be described 

by 



X = Q T (0- 1 ) 2 Q, 



(A2) 



where O is a diagonal matrix containing the semi-major 
and semi- minor axes along the diagonal, O -1 the inverse 
of O and Q T the transpose of Q. 

We define the PW contact function for two ellipses A 
and B as [H 



f AB (\) = \(l-\)v T AB Y- 



?AB, 



(A3) 



where tab = 
and Y = AX; 



-1 i 



— ya with ta the position of ellipse A, 
(1 — \)X A ~ 1 - Perram and Wertheim (58) 
has proven that /ab(A) is strictly concave on the interval 
[0, 1] and has a unique maximum Fab at A = A. Then 
Fab has the following properties: 



Fab > 1 if A and B do not overlap, 

Fab = 1 if A and B are externally tangent, 

Fab < 1 if A and B overlap. 



(A4) 



Note that Jab (A) is indeed a rational function of A and 
hence its maximum can be readily found using a Newton- 
Raphson method [62]. When Fab and A are known, we 
can scale the two ellipses by a common factor j47j so 
that they are externally tangent. The contact point is 
then determined by 



rc = r A + (1 - A)X} 



(A5) 



where n = Y _1 r^B is the unnormalized common normal 
vector at the contact point. Then re and n are used to 
compute the first-order time derivative of Fab , which has 
been derived in Ref. [47} and reads 



Fab — 2A(1 — A)n T vc, 



(A6) 



where vc is the projection of the relative velocity at the 
contact point 

vc = (v# + oj b v B c) ~ (va + u A vac), ( a ?) 

with v and uo the translational and angular velocities, 
and 



(A8) 



In principle, the contact time of two moving ellipses 
can be then numerically determined by using a Newton- 
Raphson method (62[. However, some numerical prob- 
lems are probably met and some tricks must be employed. 
For these details, the reader is referred to Ref. |46l l47j. 
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FIG. 18: For hard ellipses with k — 2 and 4, effect of the 
system size N on (a) reduced pressure, (b) nematic order 
parameter, (c) translational diffusion constant and (d) rota- 
tional diffusion constant. In (a), the results for k = 4 is shifted 
up by 5 for clarity, and the error bars in (b) are standard er- 
rors of the sample means. The system shows an isotropic 
phase in the whole density range in the case of k = 2, while 
the I-N transition occurs at 6 ~ 0.7 for k = 4. 



Once the contact time has been determined, one 
can follow the procedure of a standard EDMD simula- 
tion 63]. To accelerate the simulation, we use a neigh- 
bor list (NL) method. Note that the NL method for 
hard ellipses is more complicated than that for hard 
spheres [46|, |47|, but the essence remains unchanged. 



Appendix B: Effect of the system size 

In this section, we assess the effect of the system size. 
As introduced in Sec. I, an analysis of the system-size 
dependence of the orientational order parameter is very 
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FIG. 19: Effect of the system size N on the angular correlation 
function #2(7) for k = 4 at <j) — 0.7 and 0.8. 
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useful especially in the 2D nematic phase. However, we 
also want to explore the influence of the system size on 
other quantities in order to confirm the validity of the 
conclusions given in this work. To this end, we also per- 
formed simulations of hard ellipses with N = 200 and 
800 for k = 2 and 4. We have concluded in Sec. II that 
the system shows an isotropic phase at area fractions 
up to (j) = 0.8 in the case of k = 2 during the com- 
pression run, while the isotropic liquid transforms into 
the nematic liquid crystal at <p ~ 0.7 for k = 4. We 
find from Fig. 18 that the system size does not affect 
the above conclusion. In addition, all properties, includ- 
ing the reduced pressure the nematic order parameter 
pmax anc j diffusion constants D, don't depend on 
the system size within the statistical accuracy at least 
for the studied particle numbers. In particular, P™ ax is 



also independent of N even in the nematic phase of hard 
ellipses. This immediately indicates that the computed 
pmax j iere ^ different from the order parameter q ob- 
tained from the tensor order parameter, since previous 
simulations do show that q decrease algebraically with 
TV [23|. Thus, we cannot determine the point of the I- 
N transition from the dependence of the system size on 
pmax ^ although we have shown that P™ ax is more sensi- 
tive to the onset of the orientational order than q at fixed 
particle number. Moreover, as shown in Fig. 19, the an- 
gular correlation function #2 (r) decays faster for larger TV 
in the nematic phase, but this feature seems to become 
less evident and even disappears as (j) gets larger. How- 
ever, the results confirm again that the nematic phase in 
hard ellipses does have only quasi-LRO. 
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